Oscillatory Behavior of the Rate of Escape through an Unstable Limit Cycle 
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Suppose a two-dimensional dynamical system has a stable 
attractor that is surrounded by an unstable limit cycle. If the 
system is additively perturbed by white noise, the rate of 
escape through the limit cycle will fall off exponentially as 
the noise strength tends to zero. By analysing the associated 
Fokker-Planck equation we show that in general, the weak- 
noise escape rate is non-Arrhenius: it includes a factor that is 
periodic in the logarithm of the noise strength. The presence 
of this slowly oscillating factor is due to the nonequilibrium 
potential of the system being nondifferentiable at the limit 
cycle. We point out the implications for the weak-noise limit 
of stochastic resonance models. 

PACS numbers: 02.50.-r, 05.40. +j 



A particularly interesting phenomenon is the occur- 
rence of noise-induced transitions between attractors of 
a dynamical system. Such transitions occur in chemi- 
cal physics, where the transition is a motion across a 
transition-state surface from a reactant region to a prod- 
uct region. They also occur in statistical physics, and in 
other fields where stochastic modelling plays a role EL|| . 

If the noise is white, or has a short correlation time and 
may be approximated as white, then the probability den- 
sity of the system will satisfy a Fokker-Planck equation. 
This equation governs the way in which noise-induced 
transitions occur. By the 'rate' at which a specified tran- 
sition takes place we shall mean the probability that it oc- 
curs, per unit time. At least in finite-dimensional sys- 
tems, any such rate should fall off exponentially as e, the 
noise strength, tends to zero. (In thermal applications 
e would be proportional to kT .) In fact, each transition 
should be characterized by an activation energy AW, 
with the transition rate falling off to leading order as 
e -AW/e Computing the pre-exponential factor requires 
a careful analysis of the Fokker-Planck equation |2p|| . 

Most work has focused on the case when the compet- 
ing attractors of the dynamical system are separated by a 
separatrix (i.e., a 'ridge') containing a saddle point. How- 
ever, models where the separatrix is instead an unstable 
limit cycle arise in the context of chemical reactions con- 
strained to occur far from equilibrium Jj). Also, tran- 
sitions across an unstable limit cycle separating steady 
states of periodic vibration occur in models of stochas- 
tic resonance in bistable continuous systems [QJ|] . A full 
analysis of noise-driven escape through an unstable limit 
cycle accordingly seems called for. 

Previous work on noise-driven transitions in models 
with an unstable limit cycle is found in Refs. P,p[dl3[. 
Graham and Tel |j 10|, and Day |i~l]-|l3|], have noted that 



the nonequilibrium potential W, as a function on the 
state space of the system being modelled, will be non- 
differentiable ('wild') near the limit cycle if the system 
fails to satisfy a form of detailed balance. Naeh et al. g] 
began an analysis of the Fokker-Planck equations associ- 
ated to such models, by a method of matched asymptotic 
expansions, but their analysis assumed that W was dif- 
ferentiable at the limit cycle. 

In this Letter we begin an asymptotic analysis of the 
rate of escape through an unstable limit cycle that in- 
corporates the insights of Graham and Tel, and of Day, 
and obtain a striking result. We show that generically, 
in two-dimensional models with an unstable limit cycle 
enclosing an attractor, the rate of escape R is given by a 
non-Arrhenius formula of the form 



R~ const x e b G(|loge|)e- AW7e 



(1) 



in the weak-noise (e — » 0) limit. Here b is model- 
dependent, and the factor G(|loge|) is a model- 
dependent periodic function of | log e | . The presence 
of such slowly oscillating factors in the expressions for 
noise-dependent transition rates has not previously been 
suspected. It indicates that even in bistable dynamical 
systems with effective dimensionality as low as two, re- 
laxation phenomena may be more complicated than is 
commonly believed. Our analysis applies whenever the 
system is truly two-dimensional, i.e., is nonseparable. 

Models. — We consider models with dynamics that are 
those of a Brownian particle moving in a drift field, i.e., 
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(2) 



Here x = (a; 1 , a; 2 ) is a pair of state variables, and the 
drift field u = (u 1 ,^ 2 ) specifies the dynamics in the ab- 
sence of noise. (771, 772) is a pair of white noise processes, 
satisfying (r) a (s)r)0(t)) = 5 a p5(s — t). cr = (o % a ) is a so- 
called noise matrix that is allowed to be state-dependent 
(a 'zweibein' field). The probability density p = p(x,t) 
of such a system satisfies the Fokker-Planck equation 

p = -O = (6/2)9,9, [W (x)p\ - d t [u l (x)p] , (3) 

where the diffusion tensor D = (D 13 ) = crcr 1 . The opera- 
tor C* is the (forward) Fokker-Planck operator. We con- 
sider here the case when there is a point attractor S in 
the (x , a; 2 )-plane, with domain of attraction Q, for which 
the boundary d£l is an unstable limit cycle. 

This framework is sufficiently general that it can ac- 
comodate two-dimensional models with overdamped dy- 
namics, or one-dimensional models with underdamped 
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FIG. 1. The unstable limit cycle dQ, of the van der Pol 
model, and the MPEP, which emerges from the attrae- 
tor (0, 0) and spirals into it. The trajectories exiting from Q 
are optimal trajectories that are perturbations of the MPEP. 



dynamics. In the latter case one of the state variables 
(x , say) would be a position, and the other a velocity. 
Our simulations below are of a model of this sort, namely 



x = v 

v = -x + (v 2 - l)v + € 1/2 rj 2 {t). 



(4) 
(5) 



This is a time-reversed van der Pol oscillator, the asymp- 
totic analysis of which was begun by Day Here 
(a; 1 ,^ 2 ) = (x,v), S = (0,0), and D = diag(0, 1) is de- 
generate. The unstable limit cycle is shown in Fig. [j]. 

The Analysis. — To estimate the rate of escape 
through dQ, we use the Kramers flux-over-the-barrier 
technique [T^ ]. Suppose that escaping Brownian par- 
ticles are re-injected at the attractor S, and a steady 
state has been reached. The probability density in this 
state, which we denote po, will satisfy £*po = 0. When 
the noise strength e is small, po will be tightly peaked 
near S. po{x), at points x on the limit cycle dQ and 
outside it, will be suppressed by a factor ~ e -AW/e Te \ a _ 
tive to po(S). Since the Fokker-Planck equation has the 
form of a continuity equation, with current density J 1 [p] 
equalling pu % — (e/2)dj \_D % i p] , the escape rate R may be 
computed as the flux of probability through dQ, i.e., 



R 



mi 



J[p ] -nd£ 



po d 2 x. 



(6) 



Here n denotes the outward normal on dQ. 

To derive the oscillating formula of eq. (Q) from eq. (JTJ) , 
we introduce a WKB approximation to the steady-state 
density po in the weak- noise limit J^do)- We write 

po(x)~K(x)exp(-W(x)/e). (7) 

Here W(x) is an 'activation energy' controlling noise- 
induced fluctuations from the attractor to the vicinity 
of x. Though ([?]) resembles a Maxwell-Boltzmann distri- 
bution, W is a nonequilibrium potential, since the steady 
state is not necessarily an equilibrium state, in that 
it does not necessarily satisfy detailed balance. We set 
W(S) = 0, so that AW, the falloff rate of the escape rate, 
is the value of W attained on the unstable limit cycle. 



Substituting (fjj) into £*po = and separating out the 
0(e~ 1 ) terms yields the eikonal equation 



where 



H{x\dW/dx l ) = 0, 



H{x\ Pi ) = -D lj {x)p l p 1 +u l {x)p i 



(8) 



(9) 



is a so-called Wentzell-Freidlin Hamiltonian |Tj|. Equa- 
tion (||) has the form of a Hamilton-Jacobi equation, with 
W a classical action at zero energy. To compute W(x) 
one may simply use Hamilton's equations of motion to 
generate the zero-energy classical trajectory from S to x. 
The quantity W(x) will necessarily equal J p-dx, the line 
integral being taken from S to x along the trajectory. 
We stress that p = 'VW here is not a physical momen- 
tum: it is a mathematical artifact. Hamilton's equation 
x 1 = D^pj + u % reveals that p measures the extent to 
which the classical trajectories move against the drift it. 
Deterministic (no-noise) trajectories have p = 0. 

By separating out the O(e ) terms in £*po = one can 
show that the pre-exponential factor K(x) satisfies [0,0 



K = - (V • u + K, 



(10) 



the time derivative being a derivative with respect 
to transit time along the zero-energy classical trajec- 
tory. Here W%j = dWjdx % dxK By differentiating the 
Hamilton-Jacobi equation one can show that the matrix 
(Wjj) satisfies a Riccati equation along the trajectory: 

w,ij = /r'n ,v,u /, - u\ t w M - u k ^w M - u\ ijPl . 

(11) 

This facilitates the computation of K. 

The zero-energy classical trajectories emanating from 
the attractor, sometimes called optimal trajectories, have 
a direct physical interpretation: they are the most prob- 
able fluctuational trajectories. If a noise-induced fluc- 
tuation from S to x occurs, in the limit of weak noise 
it should occur with increasing likelihood along an opti- 
mal trajectory terminating at x. Such trajectories have 
been seen experimentally pq| . In the weak-noise limit 
the most probable escape path (MPEP) will be the least- 
action optimal trajectory extending from S to the limit 
cycle dVL. Normally this trajectory will spiral into dQ, 
rather than crossing dQ in finite time, for the follow- 
ing reason. If the MPEP crossed the limit cycle, the 
crossing point would be a 'hot spot' through which es- 
cape would preferentially occur. The tangential deriva- 
tive d t W {i.e., the tangential momentum p t ) would nec- 
essarily be zero there. But the normal drift velocity u n 
equals zero on dQ. So the second term in eq. (S) would 
vanish at the hotspot. If D is nondegenerate, eqs. (||)-(|]) 
imply that p — there, i.e., x = it. That is, the osten- 
sible MPEP would be a deterministic trajectory, which 
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is impossible. The MPEP normally spirals into the un- 
stable limit cycle even when D is degenerate. In Fig. |l| 
we show the MPEP of the van der Pol model (@)-(|). 

Multivaluedness . — Dykman, Millonas, and Smelyan- 
skiy |l7j and the present authors |4| have stressed that 
W and K may be multivalued functions of the system 
state x, since any given point x may be the endpoint of 
more than one optimal trajectory. This normally hap- 
pens in models with an unstable limit cycle, as Fig. |l| 
shows. Optimal trajectories that are perturbations of the 
MPEP do not spiral into the limit cycle. Rather, they 
approach it, and wind around the region fl a number of 
times, all the while deviating farther from the MPEP. 
They eventually exit from fl (if the perturbation is in 
the outward direction) or move back toward the attrac- 
tor (if the perturbation is inward). As a consequence, 
any point x near the unstable limit cycle is the endpoint 
of any of an infinite, discrete set of optimal trajectories, 
which differ from each other in their winding number I, 
which may be arbitrarily large. W and K are infinite- 
valued, and the WKB approximation (Q) generalizes to 
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(12) 



In the weak-noise limit, this sum is dominated by the 
term with minimum W^{x). Equivalently, fluctuations 
to any point x in fl, in the limit of weak noise, proceed 
preferentially along the physical optimal trajectory from 
the attractor to x: the least-action one. Subdominant 
trajectories contribute at larger noise strengths, however. 

We shall use (|l^) to compute the flux of Brownian 
particles over the barrier dfl. We first approximate 
and near dfl by extending the results of Nach 
et al. 0. As a first attempt, suppose that W is single- 
valued and (to leading order) quadratic near dfl, so that 
it can be approximated as AW + W. nn n 2 /2. Here n is 
the distance in from dfl, in the normal direction, and the 
second normal derivative W inn = dp n /dn < depends 
on position along dfl. The matrix equation (O) yields 



W, 



-D nn (W, r 



2u n . n W, 



(13) 



a Riccati equation along dfl. The final term in ([n]) 
has dropped out, as p = 'VW is zero on dfl if W be- 
haves quadratically there. W tTln as a function of position 
along dfl may be computed from (13) by integration [Q. 

A problem with this approach was pointed out by Gra- 
ham and Tel |^|nj. Assuming that W is single- valued 
near the unstable limit cycle is much the same as assum- 
ing that optimal trajectories that are perturbations of 
the MPEP, as well as the MPEP itself, spiral into the 
limit cycle. What actually happens near dfl is revealed 
by a Poincare section. Suppose we choose some point 
on dfl, and plot the pair (n,p n ), i.e., normal displace- 
ment and normal momentum, for each optimal trajec- 
tory that passes nearby. If W were quadratic in n, i.e., 



(0,0) 



Pn=W, nn n 



FIG. 2. A Poincare section. This sketch shows the points 
(n,p n ) generated by the optimal trajectories passing by some 
specified point on dfl. The dots are generated by the MPEP, 
spiralling into dfl. Cf. Figs. 1-3 of Graham and Tel }!□]. 



p n = dW/dn were linear in n, the points (n,p n ) would lie 
on a line with slope W ynn passing through (0,0). What 
happens instead is shown in Fig. ||. The MPEP generates 
points that tend to (0, 0) geometrically, and lie along the 
dashed linep n = W inn n. But perturbations of it generate 
points that lie along the horizontal solid lines. 

Figure ^ can be interpreted in terms of a 'return 
map' that updates (n,p n ) whenever an optimal trajec- 
tory winds once around fl. This map will have (0, 0) 
as fixed point. For the MPEP to spiral into dfl and 
yield points along the ideal line, the linearized return 
map at (0, 0) must have (1, W.nn) as an eigenvector, with 
eigenvalue less than 1. And since deterministic (p = 0) 
trajectories that 'peel off from dfl do so geometrically, 
(1,0) must also be an eigenvector, with eigenvalue greater 
than 1. By Liouville's Theorem these eigenvalues must 
be reciprocals, so we denote them c _1 and c. With each 
turn, the MPEP decreases its distance from dfl by a fac- 
tor c, and deterministic trajectories that diverge from dfl 
increase their distance from it by a factor c. 

Figure || can now be explained. Suppose the MPEP 
intersects the (n,p n ) plane at a(l,W^ nn ). Optimal tra- 
jectories that are small perturbations of the MPEP will 
intersect it at a(l,W tnn ) + Xv, where A is the pertur- 
bation strength and v is model-dependent. We write 
v = a s (l, W. nn ) + ot u {l, 0), where a s , a u ^ in general. 
By iterating the return map, we see that after winding I 
more times, the trajectories intersect the (n,p n ) plane at 



ac 



( (1,^„„) + Aa u c'(l,0). 



(14) 



The a s term has been dropped here, since it becomes 
negligible with respect to the a u term as I — » oo. It is 
the second term in (14) that gives rise to the horizontal 
solid lines of Fig. ||, as A is varied away from zero. 

On the MPEP, the noncquilibrium potential W be- 
haves quadratically near dfl. In particular, at n = aer 1 , 
W w AW + W.nn{ac- l f/2. It follows that the Z'th 
value W( L \n) of the infinite- valued function W(n), which 
arises from trajectories that wind I times around dfl, is 



W {l \n) w W{ac- 1 ) + {dW/dn) ( n — ac 



(15) 
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= AW + W. nn {ac- l ) 2 /2 + ac- l W, nn (n - ac~ l ) 
= AW + W, nn (ac- l n - a 2 C - 2l /2) . 



Each WW, as a function of the normal distance n, is to 
leading order linear, not quadratic. This has been no- 
ticed by Graham and Tel |} 10 1, who also noted that if 
one plots the physical (i.e., minimum) value of W(n), 
one obtains a piecewise linear approximation to the ideal 
parabola AW + JU n „n 2 /2. The physical W is nondiffer- 
entiable at a sequence of points converging to n = 0. 

Oscillatory Asymptotics. — To apply the Kramers 
method, we need the prefactors KW, as well as W^ l \ 
At any x, is computed by integrating eq. ( |i~(j| ) along 
an optimal trajectory that winds I times around f2, and 
terminates at x. We must distinguish here between the 
'ideal' W jnn , which is a mathematical abstraction (the 
periodic solution of the Riccati equation (p"3|)), and the 
actual second derivatives d 2 W^ /dx" l dx 3 ' . It is the lat- 
ter that appear in @. In both tffify and Fig. |, which 
were computed on the basis of the linearized return map, 
d 2 W^/dn 2 = 0, and hence d 2 W^/dx i dx j = 0, for ev- 
ery I. Keeping higher-order terms would keep the second 
derivatives 9 2 W«/8i'W from being identically zero, 
but they would still fall to zero as dil is approached. 

It follows that when computing K near dfl, we may 
replace (|l(]) by K rj — (V • u) K. In the limit of large 
winding number I, which involves integration along a 
trajectory that spirals ever closer to dft, this yields 
K(l+l)/K(l) exp[-f(V • u)dt], the integral being 
taken once around dfl. By examination, this limiting 
quotient equals c _1 . We shall write K^ l > ~ Ac~ L , where 
A is a function of position along dfl. This n-independent 
approximation is increasingly accurate as n — > 0. 

Substituting the approximations for and 
into ( )i~2| ) yields an approximation to the steady-state 
probability density near the unstable limit cycle, i.e., 



Here q = (1 + 21ogc) /2, and the summation is periodic 
in h with period unity. Equivalently, the summation is 
periodic in log e| with period 21ogc. Substituting this 
flux density into (||), which involves an integral over dfl, 
yields R ~ const x e b e~ Aw ^ e G (|log e|), where the expo- 
nent b = q — 1, and the function G(») must have period 
2 log c. This is the promised oscillatory rate formula. In- 
terestingly, b varies continuously as the model is changed. 

Discussion. — The slowly oscillating factor G ( | log e | ) is 
related to a phenomenon discussed elsewhere ||. If the 
noise strength e is small, escape across a quadratic bar- 
rier follows the formally most probable escape path (the 
MPEP) only until it gets within an 0(e 1 ^ 2 ) distance of 
the barrier. Thereafter escape occurs diffusively, rather 
than ballistically. Since the MPEP in the models con- 
sidered here spirals geometrically into the barrier d£l, 
the point at which it gets within an 0(e 1 ^ 2 ) distance 
will cycle around Q as e — > 0, periodically in log e | . 
(Cf. Day UHH].) In fact, the period will be 21ogc. 
If the effective diffusivity varies with position along dft, 
one would expect the escape rate R to be periodically 
modulated. That is what we have shown to occur. 

We expect the phenomenon of slow oscillations is rel- 
evant to stochastic resonance in multistable continuous 
systems. Hu Gang et al. [Q have recently considered 
such systems, with the addition of time-periodic forcing 
and external noise. Steady states are then periodic at- 
tractors, separated by unstable limit cycles. In the weak- 
noise limit, the rate of noise-induced transitions should 
therefore include an oscillatory factor. 
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po (n) ~ Ae- W / £ ^c-' exp{ -W, nn [ 



ac n 



a 2 c 



72]A} 



(Recall that W :1ln < 0.) It follows that J[po] ■ n, the nor- 
mal component of the probability flux density at n = 
(i.e., through dil), is to leading order 



const x e 
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The sum over winding number I may be approximated by 
a discrete analogue of Laplace's method. As e — > 0, the 
dominant terms in the sum have I « |log e| /2 logc. Let I* 
be the greatest integer less than or equal to |log e| /2 log c, 
and let h = |log e /2 logc — I*. Changing the summation 
variable to k = I — I* allows one to approximate J[po] ■ n 
in the e — > limit, up to a constant factor, by 



-AW/e 



k— — o 



c - 2 ( fc -") exp 



(-a 2 \W, m 
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